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Valeria Ferrari,' Leonardo Gualtieri,' and Luciano Rezzolla^'^ '' 

' Dipartimento di Fisica, Universita di Roma "La Sapienza" and IN FN Rome, Italy 
^ Max-Planck-Institut fUr Gravitationsphysik, Albert-Einstein-Institut, Potsdam Germany 
^SISSA, International School for Advanced Studies and INFN Trieste, Italy 
"^Department of Physics, Louisiana State University, Baton Rouge, Lousiana USA 

We present a new method for the calculation of black hole perturbations induced by extended sources in 
which the solution of the nonlinear hydrodynamics equations is coupled to a perturbative method based on 
Regge-Wheeler/Zerilli and Bardeen-Press-Teukolsky equations when these are solved in the frequency domain. 
In contrast to alternative methods in the time domain which may be unstable for rotating black-hole spacetimes, 
this approach is expected to be stable as long as an accurate evolution of the matter sources is possible. Hence, 
it could be used under generic conditions and also with sources coming from three-dimensional numerical 
relativity codes. As an application of this method we compute the gravitational radiation from an oscillating 
high-density torus orbiting around a Schwarzschild black hole and show that our method is remarkably accurate, 
capturing both the basic quadrupolar emission of the torus and the excited emission of the black hole. 
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I. INTRODUCTION 

It is well-known that matter falling into a black hole, or or- 
biting around it, is a source of gravitational radiation. The de- 
tection of these waves, and of the imprint they carry on the ex- 
citation of the black hole quasi normal modes (QNM), would 
provide an invaluable test of General Relativity in the strong- 
field regime. This process can be studied using the theory 
of black hole perturbations, the origin of which dates back to 
1957 when Regge and Wheeler Q| showed that the axial per- 
turbations of a Schwarzschild black hole can be studied by 
solving a radial equation in the frequency domain, the angu- 
lar part of the perturbation being known in terms of a suitable 
expansion in tensor spherical harmonics. 

Among the many dedicated to this subject, two papers can 
be considered as milestones in the theory of black hole pertur- 
bations. The first one is by F. Zerilli 1 2], who derived the radial 
equation governing the polar perturbations of a Schwarzschild 
black hole, thus completing the work of Regge and Wheeler 
(see |3] for a recent review). The second one is by S. Teukol- 
sky 1 4], who achieved the formidable task of reducing the 
equations for the perturbations of a Kerr black hole to a master 
radial equation. Both the Zerilli and the Teukolsky equations 
were derived in the frequency domain. 

The theory of black hole perturbations has been developed 
largely to study the stability properties of black holes and 
to determine the gravitational signals emitted by point par- 
ticles moving near black holes on different orbits (see for in- 
stance 1 5] and references therein). In addition, it has been 
used to evaluate back-reaction effects on the orbital motion of 
a point particle moving around a nonrotating black hole 1 6 ] . 

Although it clearly represents only a toy-model, the point- 
particle approximation has been very useful to study a num- 
ber of phenomena that occur in the neighborhood of black 
holes and it has also been applied to study stellar perturba- 
tions However, if one wants to investigate realistic 
astrophysical processes, in which finite-size effects and in- 
ternal dynamics need to be properly taken into account, the 
point-particle approximation is certainly too crude. Therefore, 



several attempts have been made towards a more realistic de- 
scription. In ref. f^, for instance, the point particle orbiting 
a black hole was treated as having a spin angular momentum, 
while in ref. iToll the infalling matter was assumed to be a 
finite size sphere of dust. Furthermore, in ref. flT'l a semi- 
analytic method was developed to study the radial infall of a 
shell of dust; while this latter approach has the advantage of 
being easily generalizable to sources with arbitrary shapes, it 
has the serious drawback that it can be applied only to sources 
in which the matter configurations are freely falling onto the 
black hole. It is worth stressing that in all cases mentioned 
above the black hole is assumed to be nonrotating. 

Alternative and more general approaches to treat extended 
sources have recently been developed by Papadopulous and 
Font fl^ and by Nagar et al. fisll . by combining perturbative 
approaches with fully multidimensional hydrodynamical cal- 
culations. In ref. U^, in particular, the polar perturbations of 
a Schwarzschild black hole induced by an oscillating, thick 
accretion disk (a torus) have been studied by solving the in- 
homogenous Zerilli equation in the time domain. The source 
terms were computed self-consistently by using the results of 
two-dimensional (2D) simulations in which the relativistic hy- 
drodynamical equations were solved for a torus oscillating in 
the Schwarzschild, unperturbed background. The purpose of 
that work was to compare the gravitational signal found by 
solving the Zerilli equation, with the one found using the stan- 
dard quadrupole formalism applied to the oscillating torus. 

In this paper we develop a new approach to study the per- 
turbations of a Schwarzschild black hole excited by extended 
matter sources, by solving the inhomogeneous Bardeen-Press- 
Teukolsky (BPT) equation [14, 151 in t\\s frequency domain. 
Also in this case the source for the perturbative equations is 
the result of independent 2D hydrodynamical simulations in 
the time domain. Because of the combined use of both time 
and frequency domains approaches, we refer to this as to a 
hybrid approach. 

There are at least three good reasons why our hybrid ap- 
proach may be preferrable over a pure time-domain approach 
as the one used in L13i] . Firstly, the Regge- Wheeler and the 
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Zerilli equations (hereafter RWZ equations) cannot be gener- 
alized to rotating black holes, whereas the BPT equation ad- 
mits such a generalization, i.e., the Teukolsky equation (We 
recall that the Teukolsky equation is not separable in the time 
domain because the eigenvalues of the angular functions are 
frequency dependent; this is not the case in the frequency do- 
main where the equation is indeed separable.)- Secondly, al- 
though the BPT equation is intrinsically unstable when inte- 
grated in the time domain, its integration does not offer any 
numerical difficulty when performed in the frequency domain. 
Thirdly, it is a common experience that by removing the er- 
ror in the time integration of a system oscillating around its 
equilibrium state, a frequency-domain approach improves the 
overall accuracy. The work presented here serves therefore to 
develop a new procedure for the solution of the BPT equa- 
tion with extended matter sources, whose generalization will 
allow to reach a long-standing goal: i.e., the study of the 
gravitational-wave emission from an extended source perturb- 
ing a rotating black hole. 

While the mathematical apparatus behind our hybrid ap- 
proach is well-know, our numerical implementation is more 
innovative and to test it against an astrophysically real- 
istic and computationally complex configuration we have 
here also considered the gravitational-wave emission from a 
high-density torus oscillating in its orbital motion around a 
Schwarzschild black hole. This is the same source considered 
in ref. [13] and thus offers the possibility of comparing for the 
same test case our hybrid approach with one which is instead 
in the time-domain only. 

To this scope, we recall that the possibility that extremely 
dense and massive accretion tori form around a black hole has 
been considered in recent years, initially as a possible expla- 
nation of 7-ray bursts 1 16]. In these systems, a torus of mass 
Mt ^ 0.1 — lAf© and mean density as high as ~ 10^^ — 10^^ 
g/cm"^, forms around a black hole of a few solar masses. The 
torus can be located within a few tens of km from the hori- 
zon and could also be subject to a dynamical instability which 
would induce its accretion onto the black hole on a dynam- 
ical timescale (see |17, 18] and references therein). High- 
density tori could be generated in many astrophysical scenar- 
ios (see 1 16, 18] for a general overview): in the coalescence of 
black hole-neutron star or neutron star binaries, as shown by 
several numerical simulations |19] in the collapse to a black 
hole of a rapidly rotating supramassive neutron star L2CL .2lll . 
as a byproduct of collapsars, that is, massive stars in which 
iron core-collapse does not produce an explosion, but forms 
a black hole |22]. As discussed in 1 18], these systems can be 
created with event rates comparable to that of core collapse 
supernovae. In addition, because the tori have an intrinsically 
large quadrupole moment, even small-amplitude oscillations 
would produce gravitational signals that may be detectable 
by ground based interferometers Virgo and LIGO |18]. As 
a result, the calculation of the gravitational waves from these 
tori does not only represent a useful testbed for our novel ap- 
proach, but it also offers the opportunity for a more accurate 
calculation of the emission from a realistic and intense source 
of gravitational radiation. 

The gravitational radiation emitted by a torus oscillating 



around a black hole was first computed in ref. flS*] within the 
Newtonian quadrupole formalism. However, the quadrupole 
formalism may not be able to describe a system with such 
a high energy density, and located so close to a black hole 
horizon with sufficient accuracy. In particular, the possible 
excitation of the black hole QNM cannot be revealed by the 
quadrupole approach. The same system has been studied 
in 1 13], where the perturbed equations have been integrated 
in the time domain and no excitation of the black hole quasi- 
normal modes (QNMs) was found. Conversely, as we will 
later discuss, our hybrid method applied to the same system 
allows to catch the excitation of the black hole quasi normal 
modes, though the peak they produce in the energy spectrum 
is much smaller than that corresponding to the torus oscilla- 
tion. 

The plan of the paper is as follows: in Section HI] we re- 
view the main equations which describe the perturbations of 
a Schwarzschild black hole in the RWZ and in the BPT ap- 
proaches when these are expressed in the frequency domain. 
Furthermore, we also recall the expression of the gravitational 
energy flux computed within the Newtonian quadrupole for- 
malism. In Section|ni]we present the analytical and numerical 
setup of our approach for the solution of both the perturbative 
and general relativistic hydrodynamics equations. Finally, in 
Section llVl we discuss the results of the numerical integrations 
for the representative case of an oscillating torus, while in Sec- 
tion |V] we draw our conclusions. A number of Appendices 
contain details on the form of the source functions employed 
for the solution of the perturbative equations. 

II. A BRIEF REVIEW OF BLACK HOLE 
PERTURBATIONS IN THE FREQUENCY DOMAIN 

In what follows we briefly summarize the relevant equa- 
tions needed to describe the perturbations of a Schwarzschild 
spacetime in the RWZ and in the BPT formulations, as they 
will later be used in our hybrid approach. We will also write 
the expression of the energy flux computed within the Newto- 
nian quadrupole approximation since it will be used for com- 
parison with the corresponding quantities computed within 
the perturbative approaches. 

A. The RWZ approach 

In the RWZ approach (ll|3], the equations describing the ra- 
dial behaviour of the perturbations of a Schwarzschild black 
hole are reduced to two wave equations after expanding all 
perturbed tensors in tensorial spherical harmonics and after 
separating the radial from the angular part. In particular, they 
are expressed for two suitable combinations of the compo- 
nents of the perturbed metric tensor and which are referred 
to as the Zerilli function, Z^'^^, and the Regge- Wheeler func- 
tion, respectively. Hereafter we will indicate with a 
"tilde" all the perturbative variables that have a time depen- 
dence which is not explicitly indicated and with i, m the in- 
dices of the spherical harmonic decomposition. 
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Since the system we will consider, namely the oscillating 
torus perturbing the spherically symmetric black hole space- 
time, is axially symmetric, the perturbations are also axi- 
ally symmetric, and hence in what follows we will consider 
m = 0, with the Zerilli and Regge- Wheeler functions ex- 



pressed as (r, i), ' (r, t), with l>2. 
Their Fourier transforms are simply given by 
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Z)^'{r, t)dt , 



and satisfy the inhomogeneous wave equations 



(2.1) 



(2.2) 



where the upper index "+" refers to the Zerilli equation, 
the "— " one to the Regge-Wheeler equation, and the symbol 
IL'-*'' is just a shorthand for the differential operator 



1L(±) EE dl 



V''-- 



Note that in order to provide mathematically consistent 
boundary conditions at the event horizon, the latter is moved 
to negative spatial infinity in terms of the "tortoise" radial co- 
ordinate = r + 2M In (r/2A/ — 1), and that the Zerilli and 
Regge-Wheeler potentials in equation M.2\ have explicit ex- 
pressions given respectively by 



2(r - 2M) 



r4(Ar + 3M): 



■[A2(A + l)r3 



-3Af AV^ + OAf^Ar + OAT^] 



_ 2(r- 2A/) 



[(A + l)r - 3M] 



(2.4) 
(2.5) 



Here A = - 1)(^ + 2)/2, while S^^\r,uj) are suitable 
combinations of the components of the stress-energy tensor 
of the extended matter source, after having been expanded in 
tensor spherical harmonics and Fourier transformed in time. 
Explicit expressions for the source functions S^^\r,uj) can 
be found in Appendix IA II 

Equations M.2\ can to be solved numerically by imposing 
outgoing boundary conditions at spatial infinity {i.e., r, r» — > 
oo), and ingoing boundary conditions at the event horizon 
(i.e., r = 2M, r» —oo). The energy spectrum of the 
gravitational radiation measured by a distant observer and in- 
tegrated over the solid angle, can then be expressed in terms 
of the asymptotic amplitude of zi^"^ and zi ^ as 



dE 



oo 

E 

oo 

= E 

1=2 



duj 

1 (^ + 2)! 
(£-2)! 



. (2.6) 



B. The BPT approach 

In the BPT approach (30. the perturbations of the Weyl 
scalar S'i>4, which describes the outgoing gravitational ra- 
diation, are expanded in spin-weighted spherical harmonics 



_2<5'|„j and integrated over the soUd angle dfl = sin ( 

^im{r, <) = f dQ -2S}^{d, 4>)rH^4{t, r, i?, 0) . (2.7) 



Also in this case, being the perturbing stress-energy tensor ax- 
ially symmetric the only nonvanishing components are those 
with m — Q and hence perturbed Einstein equations take the 
form of inhomogeneous wave equations 



IL ^i,{r,uj) = -S^ (r,cj), 



(2.8) 



where (r, uj) is the Fourier transform of ^^(r, t), the sym- 

BPT 

bol IL is just a shorthand for the differential operator, 
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IL EE /^^dr-^dr + U{r) 



(2.9) 



(2.3) with A = r{r— 2M) and the potential U having explicit form 



U : 



r^Lj^ + 4i(r - M)r^uj 



2A 



(2.10) 



As for the RWZ equation, the source terms ir,uj) are 
suitable combinations of components of the stress-energy ten- 
sor after having been expanded in spin- weighted spherical 
harmonics and Fourier transformed. Explicit expressions for 
the source functions S^^^ [r, uj) can be found in Appendix IA 21 
Equation ( 12. 8> can be integrated numerically imposing the 
same boundary conditions discussed for the RWZ equation. 
In this case, the energy spectrum of the gravitational radiation 
collected by a distant observer and integrated over the solid 
angle is expressed as 



dE""^^ ^ dE, 



duj 



E duj 
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e=2 
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(2.11) 



C. The Newtonian quadrupole approximation 

The Newtonian quadrupole approximation is based on the 
assumptions that the gravitational field is weak, that the veloc- 
ity of matter in the source is much smaller than the speed of 
light, and that the source itself is small compared to the char- 
acteristic wavelength of the emitted gravitational radiation. In 
this approach, then, the amplitude of the gravitational wave is 
estimated in terms of the quadrupole tensor associated to the 
source 



T°°(t,x)x'x^d^x = 1,2,3). (2.12) 



where T is the stress-energy tensor of the matter source. The 
resulting gravitational-wave energy-spectrum is then given by 
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'J2\Qvi^)Q^li^)\ (2-13) 
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where Qijiuj) is the Fourier transform of the traceless 
quadrupole tensor Qij 

Qtj = 9y - ^^ij X! '^'^^'^^ ■ (2-14) 

kl 

As we will comment more extensively in Section Hvl we 
will compare expressions i2.6\ . (I2.11> and (I2.13> for the grav- 
itational emission coming from an oscillating torus and point 
out the respective advantages and disadvantages. 

III. THE HYBRID APPROACH 

We will now give a general and brief account of our hy- 
brid approach, sketching its most relevant parts and its numer- 
ical implementation. We will first discuss how to construct 
the components of the stress-energy tensor, needed as source 
terms in the perturbative equations, through the solution of the 
2D hydrodynamic equations in a Schwarzschild background. 
We will then illustrate how to find the solution of the RWZ and 
BPT equations once the sources are known in the frequency 
domain. 



A. Computing the sources from the hydrodynamic equations 

As mentioned in the Introduction, any realistic and ex- 
tended matter source whose dynamics goes beyond the el- 
ementary free-fall from infinity, will need to be described 
through the solution of the general relativistic hydrodynamics 
equations on the background black hole spacetime. In partic- 
ular, this is accomplished by computing at the different space- 
time points covered by a numerical grid, the components of 
the stress-energy tensor T^,y{t, r, 9) necessary for the calcula- 

— ~BPT 

tion of the sources and 

For a matter source represented by a perfect fluid with four- 
velocity u and described by the stress-energy tensor 

Tf^iy = (e + p)Uf^u^ + pg^_,i, = phUf_,u^ + pgf_,^ , (3.1) 

where g^i, are the coefficients of the metric which we choose 
to be Schwarzschild. Here, e, p, p, and h — {e + p)/p are 
the proper energy density, the isotropic pressure, the rest-mass 
density, and the specific enthalpy, respectively. In practice, all 
of these quantities are computed by the solution of the conser- 
vation equations for the stress-energy tensor and baryon num- 
ber density 

V^T^"" = , (3.2) 
V^ipun = , (3.3) 

where V indicates the covariant derivative in the background 
Schwarzschild spacetime, together with an equation of state 
(EOS) relating the pressure to other thermodynamic al quanti- 
ties. For simplicity we will hereafter model the fluid as ideal 
with a polytropic p = up"' = pe(7 — 1), where e — e/ p ~ 1 
is the specific internal energy, k is the polytropic constant and 
7 is the adiabatic index. 



In order to preserve their conservative nature, we cast equa- 
tions i3.2\ . j3.3l l in the form of a flux-conservative hyper- 
bolic system after introducing suitable "conserved" variables 
rather than in terms of the ordinary fluid, or "primitive ", vari- 
ables. In this case, equations ( l3.2> -( ll!2t assume the form L23ll 

9U(w) ^ 9[7^F''(w)] ^ 
dt dr 

"'^'"■"' =S(»). <3.4) 

where U(w) = {D, Sr, Se, S^),F^ and S are the state-vector, 
the fluxes and the sources of the evolved quantities, respec- 
tively (see fvf] for the explicit expressions for F* and S in a 
Schwarzschild spacetime). The following set of equations 

D = pT, 

Sj = phT^Vj , (3.5) 

together with the ideal-fluid EOS provide the relation be- 
tween the conserved and primitive variables in the vector 
w = {p,Vi,e). Here F = au* = (1 — where 

= jijV^v^ is the Lorentz factor measured by a local static 
observer Note that the covariant components of the three- 
velocity are defined in terms of the spatial 3-metric to be 
Vi = 'jijV^ , where = /au*' . (Although in axisymmetry, 
we evolve also the azimuthal component of the equations of 
motion, so that the index j takes the values j = r,9,(p.) 

The numerical code used in our computations is 
the same used in ref. jTs^l and it performs the nu- 
merical integration of system ( 13 .4> using upwind high- 
resolution shock-capturing (HRSC) schemes based on ap- 
proximate Riemann solvers. Exploiting the flux conserva- 
tive form of equations (13. 4> . the time evolution of the dis- 
cretized data from a time-level n to the subsequent one 
n + 1 is performed according to the following scheme 

^ (fL+1/2 - ^1-1/2) + A< S.,, , (3.6) 

where the subscripts i,j refer to spatial (r, 9) grid points, so 
that U"^ = XJ{ri,9j,t'^). The inter-cell numerical fluxes, 

j ^iid j±i/2' computed using Marquina's ap- 
proximate Riemann solver il8.1 . A piecewise-linear cell re- 
construction procedure provides second-order accuracy in 
space, while the same order in time is obtained with a con- 
servative two-step second-order Runge-Kutta scheme applied 
to the above time update. 

Our computational grid consists of Nr x Ng ~ 250 x 84 
zones in the radial and angular directions, respectively, cov- 
ering a computational domain extending from r„iin = 2.1 to 
''max = 30 and from to tt. The radial grid is logarithmically 
spaced in terms of the tortoise coordinate with the maximum 
radial resolution at the innermost grid being Ar = 6 x 10^"*. 
As in ref. L18J . we use a finer angular grid in the regions that 
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are usually within the torus and a much coarser one outside. 
The boundary conditions adopted, the treatment of the vac- 
uum region outside the torus with a low density atmosphere, 
and the procedure for recovering physical variables from the 
conserved quantities D and Si are the same as those used in 
ref. fr3l . The interested reader is referred to that work for 
further details. 



B. Solution of the RWZ and BPT equations 



Before solving eq. \2.2l with ingoing wave boundary con- 
dition at the horizon and outgoing wave boundary condition 
at infinity, we determine the solution of the corresponding ho- 
mogeneous equation 



IL 



= 



(3.7) 



for an assigned value of lo, finding two particular solutions, 
zj^-* ^j^"* which satisfy, respectively, ingoing wave 
boundary conditions at the horizon 



and outgoing wave boundary conditions at infinity 

/^(i) out/ N iuir, 

ZjI ' (r, ^ oo, = e * . 
Since they are independent solutions, their Wronskian 



(3.8) 



(3.9) 



(±) in^ 



dZ, 



(±) out ' 



,(±) in 



(3.10) 

is constant. The solution of the inhomogeneous equations 
(|^}, zi^\ is then found as a convolution integral of the 



source with 



(±) 



and Z, 



(±) out 

e 



{r,uj) 



,(±) out 

'_e 

,(±) out 

'_e 



y{±) in M±) 12 



2M 



y{±) out r-(±) 12 



(3.11) 

At spatial infinity the Zerilli function can be written as l24ll 



(±) 



where 



.)2M 



y(±) in „(±) 2 

dr^ 



(3.12) 



(3.13) 



In this way, the functions A^^^^ and A\ ^ represent the ampli- 
tudes of the gravitational radiation with polar and axial parity 
and, as expressed by ( 12. 6> . they both contribute to the gravita- 
tional spectrum 



dE 



dw 



y- 



1 (^ + 2)! 



1=2 



8 {(.~2)\ 



(3.14) 



As done for the RWZ equations, also for the BPT equation 
( 12. 8> we first solve the corresponding homogeneous equation 



BPT 

IL ■^i{r,u})^Q 



(3.15) 



for an assigned value of uj and find two independent solutions 
which satisfy, respectively, the condition of a pure ingoing 
wave at the black hole horizon and of a pure outgoing wave at 
infinity 



-oo, Uj) 



(3.16) 



(3.17) 



where the correcting factors A and in ( I3.16> and ( I3.17l i 
are the result of the asymptotic expansion of the BPT equa- 
tion fl4','T?]. 

The solution ^'^(r, uj) of the inhomogeneous equation ( I2.8l l 
has a form analogous to the one for the RWZ equations ( 13.1 H 
and also in this case for — > oo it can be expressed as 



'^i{r^,uj) At {uj)r^e 



3 Aujr^ 



with 



A (^) 



1 



dr'-' 



(r,u;) 



W J2M A2 

where the Wronskian is now defined as 



-4 



\ dr 



dr 



(3.18) 



(3.19) 



(3.20) 



Also in this case, the complex amplitude A^ describes the 
gravitational-wave amplitudes with polar and axial parity, so 
that the gravitational spectrum ( 12. lU becomes 



dE 



du) 



E 

1=2 



At 



(3.21) 



It is important to note that the solutions v|/^"/°"' of the ho- 
mogeneous equations ( I3.16> and (I3.17> . can be expressed in 



terms of the solutions Z, 



( — ) in/out 

e. 



of the corresponding ho- 
mogeneous RWZ equations (13. 7> through the Chandrasekhar 
transformation ll25ll 



in/out 



'(^ + 2)! 



in/out 



r - 3M 



A 



j,2 ^tr 



{— ) in/out 



( — ) in/out 



(3.22) 



In what follows we will use this transformation to avoid the 
resolution of (I3.15>. 
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C. Combining the two approaches 

Having described the distinct approaches for the calcula- 
tion of the sources in the time domain and the solution of 
the perturbative equations in the frequency domain, we now 
discuss how to combine the two methods within our hybrid 
approach. Since the methodology applies unchanged whether 
we are considering the solution of the RWZ or of the BPT 
equations, we will drop this distinction in what follows. 

Assume therefore that the solution of the hydrodynamical 
equations i3.2\ . i3.3\ as illustrated in Sect. IIII Al has provided 
the components of the stress-energy tensor T^i/(i, r, 9) for all 
the spacetime events that are relevant. The first step for the ef- 
fective calculation of the sources consists then in the removal 
of the angular dependence. This is done by calculating for 
the different values of £, the integrals over the 2-sphere of 
the stress-energy tensor components, as given in eqs. (|^n)- 
(IA17> and (IA33> - (IA35> to compute the following quantities 

Ae, aI'\ Bf \ Be, Qe, Fe, 5,, 

T{n)(n)e, T(^n)(m)e, T(^m)(nf?i)e ■ (3.23) 

As a result, all of the above quantities are functions of {t, r) 
only and are evaluated at the discrete spacetime points r^, tj, 
where coincide with the radial gridpoints of the 2D code. 
The time levels tj — jAt e [0, T], on the other hand, are 
chosen so that T ^ t, where r is the typical timescale for 
the problem (e.g., the torus oscillation timescale for the prob- 
lem at hand and T/r ~ 100), and Ai < r (e.g., At/r ~ 
10^^ — 10^ ■^). We note that in general the solution of the hy- 
drodynamical equations does not take place on equally spaced 
spacelike hypersurfaces since the time step is automatically 
adjusted on the basis of the dynamics of the matter source. As 
a result, the evaluation of the functions J3.23> on the equally 
spaced time levels tj requires in general an interpolation pro- 
cess. 

Once the interpolated timeseries for the functions (I3.23> 
have been constructed at each gridpoint r^, these are Fourier 
transformed yielding the corresponding frequency-domain 
functions 

Ae, bI°\ Be, Qi, F,, De, 

T{n)(n)e.^ T(^n)(m)ej T(^m)(nf?i)e (3.24) 

which are evaluated at the points (r^, cjj), with ujj = jAuj = 
j27r/T. The functions (I3.24> are then suitably combined as 
in equations (IA20> . (IA21> and (IA32> to provide the com- 
pact forms for the source terms in the frequency domain. 

The next step consists of integrating the homogeneous 
equations \3.1\ to find the two independent solutions ( 13. 8> . 
( 13. 9> . The numerical solution is found using a 4*^-order 
Runge-Kutta integrator with adaptive stepsize and with 
boundary conditions given by expressions (13. 8> and (13. 9> . Us- 
ing the newly computed RWZ solutions z''^^ ™^°^^{ri,ujj) 
and the transformation (I3.22> . it is possible to also compute 
the two independent solutions of the homogeneous BPT equa- 
tion ( I3.15> . \I'^"^°"'(ri, Wj). Finally, the amplitude of the 



emitted gravitational wave can be computed by numerically 
evaluating the integrals (I3.13> and ( I3.19> . which then yield 
the energy spectra ( I3.14> and ( I3.21> . 

With the exception of the handling of the solution of the rel- 
ativistic hydrodynamical equations, the procedure described 
so far for the combination of the time and frequency-domain 
approaches is straightforward and with minimal computa- 
tional requirements. However, great attention needs to be paid 
in to avoid unphysical results. In particular, even when the 
timeseries extend over several tens of dynamical timescales 
and the time sampling is also very high, the calculation of the 
Fourier transforms can be inaccurate and this can be problem- 
atic particularly at very large frequencies. Indeed, we have 
found that the energy spectrum can become divergent at fre- 
quencies above a few kHz even when the hydrodynamical 
evolution is carried out over ~ 100 dynamical timescales. 

The reason for this is most easily seen within the Green's 
function approach; in this case, in fact, when u) increases 
above a few kHz, the Wronskian of the homogeneous solu- 
tions tends to zero very rapidly. At least in principle, this 
rapid decay should be compensated by an equivalent decrease 

i^iz) BPT 

of the source function S\ or Sg such that the wave am- 
plitude would remain finite. In practice, however, this is not 
necessarily the case and if the decay in the source functions 
is not sufficiently rapid, because for instance their values at 
high frequencies are not sufficiently accurate, this would in- 
evitably lead to the high-frequency divergences we have ob- 
served. Fortunately, a simple solution to this otherwise serious 
problem is possible. In the limit of an infinite time integration 
interval, in fact, it is possible to replace the Fourier transform 
at a given frequency lo of a. timeseries h{t) with the Fourier 
transform of its time derivative divided by uj, i.e. 

This identity, which is strictly true if h{t) — for t ±oo, 
can be exploited to compensate the loss of accuracy in the 
Fourier transform. As a result, in our approach, and before 
Fourier transforming the source terms Ag, A^^' , ... in ( I3.23t . 
we compute their first and second time derivatives, evaluate 
the corresponding Fourier transform and finally divide the re- 
sult by a>^. This method effectively removes the divergence 
appearing at large frequencies and has proven to work well 
also in the case in which the timeseries is, in practice, finite. 

IV. A REPRESENTATIVE EXAMPLE: AN OSCILLATING 
HIGH-DENSITY TORUS 

The procedure discussed in the previous two Sections is to- 
tally generic and could in principle be used to investigate the 
perturbations induced on a Schwarzschild spacetime by an ex- 
tended matter source when the latter dynamics is simulated 
consistently in more than one spatial dimension. However, 
as mentioned in the Introduction, we will here concentrate on 
a proof of its effectiveness by considering the perturbations 
induced by a non-selfgravitating torus orbiting and oscillat- 
ing around a Schwarzschild black hole. The latter will be 
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TABLE I: Main properties of the constant angular momentum toroidal neutron star models used in the numerical calculations. From left to 
right the columns report: the name of the model, the star-to-hole mass ratio Mt/M, the polytropic constant k, the specific angular momentum 
£ (normalized to M), the inner and outer radii of the toroidal neutron star rin and rout, the radial position of cusp rcusp, the radial position 
of the centre rcentre (all radii are in units of the gravitational radius rg), and the orbital period at the centre of the torus toib, expressed in 
milliseconds. The last two columns indicate the density at the centre of the torus and the average density of each model, respectively, both in 
cgs units. All of the models share the same mass for the black hole, M — 2.5Mq and adiabatic index 7 — 4/3. 



Model 


Mt/M 


K 


I 


nn 


rout 


^cusp 


^centre 


iorb 


pccntre 


(P) 






(cgs) 












(ms) 


(cgs) 


(cgs) 


(c) 


0.1 


0.96x10" 


3.8000 


4.576 


15.889 


4.576 


8.352 


1.86 


1.14x10" 


4.73x10" 


(e) 


0.1 


7.0x10" 


3.7845 


4.646 


14.367 


4.646 


8.165 


1.80 


1.61x10" 


6.43x10" 


(f) 


0.1 


l.OxlO" 


3.8022 


4.566 


16.122 


4.566 


8.378 


1.87 


1. 10x10" 


4.48x10" 



simulated with the 2D general relativistic code discussed in 

sect.mra 

In order to to construct the background initial model for the 
torus, which we subsequently perturb, we consider a perfect 
fluid described by the stress-energy tensor ii. li and in circular 
non-geodesic motion with four- velocity u" ~ (w',0, AjM*^) = 
M*(l, 0, 0, n), where fl — ri{r, 9) = u'^/u* is the coordinate 
angular velocity as observed from infinity. Enforcing the con- 
ditions of hydrostatic equilibrium and of axisymmetry simpli- 
fies the hydrodynamics equations considerably and for a non- 
selfgravitating fluid these reduce to Bernoulli-type equations 



P 



1 - ni 



(4.1) 



where i ~ r,9, W = In(ut) and / = —u^/ut is the specific 
angular momentum. 

The simplest solution to equations j4.U is the one with 
I = const., since in this case the equipotential surfaces can 
be computed directly through the metric coefficients and the 
value of the specific angular momentum. Note that at any 
point in the (r, 6) plane, the potential W can either be positive 
(indicating equipotential surfaces that are open) or negative 
(indicating equipotential surfaces that are closed). The case 
W ^ refers to that special equipotential surface which is 
closed at infinity (see refs. llTlflSll for details), local extrema 
on the equatorial plane of closed equipotential surfaces mark 
the radial positions of the cusp, rcusp, and of the "centre" of 
the torus, Vc- At these radial positions the specific angular 
momentum must be that of a Keplerian geodesic circular orbit 
which can effectively be used to calculate the position of both 
the centre and the cusp. 

Stationary solutions with constant specific angular momen- 
tum are particularly useful since in this case the angular ve- 
locity is fully determined as = Igtt /g^,^, and if a polytropic 
EOS is used, the Bernoulli equations ( 14. 1> can be integrated 
analytically to yield the rest-mass density (and pressure) dis- 
tribution inside the torus as 



7-1 

K7 



CXp(Win - T4^) - 1] 



1/(7-1) 



(4.2) 



1.25 - 



0.75 




where Win = W^(nn,7r/2). Clearly, different configurations 



100 1000 
(2tt)-1w (Hz) 



FIG. 1: Harmonic behaviour of an oscillating torus. Both panels 
refer to model (c) of TableQwhich was perturbed with a radial veloc- 
ity having 77 = 0.05. The left panel shows the time variation of the 
maximum rest-mass density normalized to its initial value for a repre- 
sentative portion of the timeseries. The right panel shows the power 
spectral density (PSD) calculated over 100 orbital timescales, with 
the eigenfrequency of the / mode indicated with a vertical dashed 
line at 227 Hz. 



can be built depending on the value chosen for /, with finite- 
extent tori resulting when < I < Znib, with /„ib = 4, 
^ms = 3-\/3/2 being the the specific angular momenta cor- 
responding to orbits that are marginally bound or marginally 
stable, respectively. 

The configurations built in this way have a large quadrupole 
moment but being stationary they do not produce any time- 
varying perturbation to the background spacetime. Because 
of this we introduce parametrized perturbations that would 
induce a small outflow through the cusp and excite a quasi- 
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100 1000 10" 100 1000 10* 

(27t)-1cj (Hz) (27t)-1cj (Hz) 



FIG. 2: Gravitational-wave energy spectrum produced by an oscil- 
lating, high-density torus orbiting around a 2.5 AIq Schwarzschild 
black hole. The radiation has been computed for the £ — 2 mode 
using eq. <2.13> within the Newtonian quadrupole approximation. 



harmonic behaviour in the hydrodynamical variables (the in- 
terested reader is referred to refs ll26H27ll28ll29ll for a detailed 
discussion of the harmonic properties of this type of oscilla- 
tions). More specifically, we modify the stationary equilib- 
rium configuration with a small radial velocity which we have 
expressed in terms of the radial inflow velocity characteris- 
ing a relativistic spherically symmetric accretion flow onto a 
Schwarzschild black hole, i.e., the Michel solution ll30ll . Us- 
ing T] to parameterize the strength of the perturbation, we have 
specified the initial radial (covariant) component of the three- 
velocity as 

Vr = ?7(«r)M.oh„l ■ (4-3) 

As a result of the introduction of the initial perturbation, the 
torus acquires a linear momentum in the radial direction push- 
ing it towards the black hole. When this happens, the pressure 
gradients become stronger to counteract the steeper gravita- 
tional potential experienced as the torus moves inward, thus 
increasing the central density and eventually pushing the torus 
back to its original position. The resulting oscillations are es- 
sentially harmonic in the fundamental /-mode of oscillation 
but other, higher-order p modes are also excited and appear at 
frequencies that are in a ratio of small integers 
This is summarised in the two panels of Fig.[2which refer to 
model (c) of Table |I] which was perturbed with a radial ve- 
locity having 77 = 0.05 and evolved for several tens of the 
orbital timescale torb — 1-86 ms. The left panel shows as a 
function of time the variation of the maximum rest-mass den- 
sity (i.e., the rest-mass density at the centre of the torus pc) 
normalized to its initial value for a representative portion of 



FIG. 3: Same as Fig.|2|but derived integrating the RWZ equations 
in the frequency domain for £ = 2 and using eq. 12. 6> . 



the timeseries. The right panel, on the other hand, shows as 
a function of the frequency = uj/{2tt), the power spectral 
density (PSD) in arbitrary units of the timeseries calculated 
over about 100 orbital timescales. Indicated with a vertical 
dashed line at D — 227 Hz is the eigenfrequency of the fun- 
damental (/) mode of oscillation and smaller peaks at integer 
and semi-integer multiples of D are also visible and are re- 
ferred to the first overtones (see ref. for a discussion of 
these modes). 

As the simulation is carried out and the hydrodynamical 
equations are solved along the lines discussed in Sect. IIII Al 
the source functions ( I3.23t are calculated at the different ra- 
dial gridpoints and stored as distinct output. Once the dy- 
namics has been followed for a sufficient time-span compris- 
ing several tens orbital timescales, the simulation in the time- 
domain is stopped and the values of source functions ( I3.23t 
throughout the simulated spacetime are read-in at the differ- 
ent radial positions. This marks the begin of the analysis in 
the frequency-domain, which first produces the source func- 
tions ( I3.24t and then proceeds to the solution of the perturba- 
tion eqs. i2.2\ and ( 12. 8> along the lines discussed in Sect. IIII Bl 
The results of the hybrid approach are then summarised in 
Figs. |2j0 which show the energy spectrum of the emitted 
gravitational radiation as a function of the frequency. 

The first figure, in particular, has been computed within the 
Newtonian quadrupole formalism through expression j2.13t 
and being a faithful mirror of the hydrodynamical quantities, 
it shows a main peak at 227 Hz as well as the smaller peaks al- 
ready discussed for the right panel of Fig. ^ Note that despite 
the use of the procedure described at the end of Sect. IIII CI 
the very steep dependence on the frequency as oj^ of the en- 
ergy spectrum [cf. eq. ( I2.13H . produces an incorrect but finite 
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100 1000 10* 

(27t)-1cj (Hz) 



FIG. 4: Same as Fig.|2|but obtained integrating the BPT equation in 
the frequency domain for i = 2 and using eq. <2.1 1> . 



growth at high frequencies, where the accuracy in the calcu- 
lation of the sources is smaller. Of course, being rooted in 
a Newtonian approximation, the energy spectrum in Fig. ^ 
shows no sign of a contribution coming from the QNMs of 
the black hole that, for a mass of i\/ = 2.5 M© are expected 
at = 4.828 kHz for the the fundamental £ = 2-QNM. 

Fig- E shows instead the equivalent gravitational-wave en- 
ergy spectrum resulting from the hybrid approach employ- 
ing the solution of the RWZ equations [cf. eq. ( 12. 6H for 
£ = 2. A number of interesting features should be noted. 
Firstly, at low frequencies the spectra derived through the per- 
turbative approach shows a rather good agreement with the 
quadrupole spectrum given in Fig.|2] This confirms what al- 
ready noted in fisll and that, at least for this type of sources, 
the quadrupole approach captures the most important qualita- 
tive features of the energy spectrum. Secondly, the spectrum 
is not diverging, nor growing at large frequencies as a result 
of the use of the transformation j3.25l l. Finally and most im- 
portantly, the spectrum is not monotonically decreasing but 
shows a distinctive and broad peak at v ^ 4.4 kHz, rather 
close to the position of the expected fundamental £ — 2 QNM 
of the black hole and which we interpret as the excitation of 
the QNMs of the black hole resulting from the perturbations 
induced by the oscillating torus. Although its contribution is 
energetically very small (i.e. less that 2 x 10^* of the to- 
tal emitted energy is produced at frequencies larger than 2 
kHz), the ability to isolate this peak is of great importance 
to assess the validity of the hybrid approach and to prove its 
effectiveness in modelling black hole spacetimes. The pres- 
ence of such a peak, in fact, is expected on the basis of simple 
considerations but was not detectable in the complementary 
work of Nagar et al. 1 13] based on a time-domain perturbative 



approach despite the corresponding PSD reached values was 
well below the ones at which our black hole peak appears {cf. 
Fig|3]with the left panel of Fig. 1 in ref. Q). 

Finally, we show in Fig. 0] the gravitational-wave energy 
spectrum resulting from the hybrid approach employing the 
solution of the BPT equation {cf. eq. ( I2.11H for ^ = 2. As 
expected, much of what found for the RWZ equations contin- 
ues to hold when the hybrid approach is computed using the 
BPT equation, namely: very good match with the quadrupole 
approximation at low frequency where most of the energy is 
emitted and the appearance of a broad peak associated with 
the black hole quasi-normal ringing. 



V. CONCLUSIONS 

We have presented a new hybrid approach to study the per- 
turbations of a Schwarzschild black hole excited by extended 
matter sources in which the solution in the time-domain of 
the relativistic hydrodynamical equations in a multidimen- 
sional spacetime is coupled to the solution in the frequency- 
domain of either the Regge-Wheeler/Zerilli equations or of 
the Bardeen-Press-Teukolsky equations. 

We believe that our hybrid approach may be preferrable 
over a pure time-domain approach based on the RWZ equa- 
tions for at least three different reasons. Firstly, the Regge- 
Wheeler and the Zerilli equations cannot be generalized to 
rotating black holes, whereas the BPT equation admits such 
a generalization. Secondly, the BPT equation is intrinsically 
unstable when integrated in the time domain but its solution 
is regular computed in the frequency domain. Thirdly, for a 
system undergoing small oscillations around an equilibrium 
state, a frequency-domain approach is expected to be signifi- 
cantly more accurate than one based in the time-domain, espe- 
cially at high frequencies, where the black hole contributions 
are expected to emerge. Within this framework, therefore, the 
work presented here serves as a first step towards the study 
of the gravitational-wave emission from an extended source 
perturbing a rotating black hole. 

As a test of its effectiveness and to provide a close compar- 
ison with what done so far in perturbative approaches based 
in the time-domain, i.e., ref. lll3ll . we have considered as ex- 
tended source to the perturbation equations the oscillations of 
a high-density torus orbiting around a Schwarzschild black 
hole. The dynamics of the torus has been followed numeri- 
cally using a 2D numerical code solving the general relativis- 
tic hydrodynamics equations using HRSC schemes in spheri- 
cal polar coordinates. 

Overall, the results of our study show that the hybrid ap- 
proach (either with the RWZ equations or with the BPT one) 
works remarkably well, producing an energy spectrum of the 
emitted gravitational radiation which contains both the ba- 
sic quadrupolar emission of the torus as well as the excited 
emission of the black hole. The latter feature was completely 
washed out in the complementary time-domain approach of 
ref. ill. 

Being totally generic, the formalism and methodology pre- 
sented here is could be used under more general conditions 
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and in particular also with sources coming from 3D numerical 
relativity codes as long as an accurate evolution of the matter 
sources is possible. Future work in this line of research will 
therefore comprise both the extension of the present results 
to 3D codes, as well as the generalization of the perturbative 
formalism to the Teukolsky equation in order to model the 
gravitational emission of rotating black holes when perturbed 
by extended sources. 
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APPENDIX A: GENERIC SOURCE TERMS 
1. Source for the RWZ equations 

The source of the RWZ equations \2.2\ . sf^'^ (r, oj), can be 
expressed in terms of the components of the stress-energy ten- 
sor T^,y{t, r, 9, (p) as follows. We first decompose T^j/ in ten- 
sor spherical harmonics for the the relevant components 



TtAr,t) 

Trr{r, t) 
Tra{r, t) 

Tabir, t) 
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im 



im 
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where a = 6,(j}, ^ab = diag(l, sin^ 9) is the metric tensor of 
the sphere, F^™ are the scalar spherical harmonics, F^™ = 

QY'^m^Q^a ^jjj ^jjjj ^ = {£ -!){£ + 2)/2, A' = e{i + 1). 
The matrices 

gim ^ ^-Y'^ /sin9,sin9Y%) , 



yim 
^ab 



Qim 

'-'ab — 



X 



im 



- sin^ 9W 



im 



A:^"/sin6l sm9W 



i9W'^ 



sin 61 a: 



im 
im 



(A6) 
(A7) 

(A8) 



are vector and tensor spherical harmonics, and collect all of 
the angular dependence. 

Next, using the orthogonality properties of 
Y^"\ Si"', S"^;", Z^^' and the axial symmetry of the 
stress-energy tensor it is possible to show that only the ra — Q 
harmonics contribute to the expansion (IA5> and that 



(All) 
(A12) 
(A13) 
(A14) 
(A15) 

. 2 r^^W^o, (A16) 
s\\r'9 J 

De = -iA22^ / d9sin9 ( ^^TeA , (A17) 



Ae^2TT J d9s\n9Y^°Trr , 
I^^^ = -i\/227r J d9 sin 9Y'^°Ttr , 
= -iAi27r / d9sm9Tt0Y^°g 



Be = Ai27r J d9 sin 9TreY'^^g , 

= Ai27r / d9 sin 9 ^Tr^Y'^\ 
I sniO 



Ff = A227r / ^6* sin 61 Tf 
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where 
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(A18) 
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In terms of the Fourier transform of these quantities, the 
source terms 5^^' (r, lo) are 

S^+\r,oj) = uqA^^^ +aiA^p + I3A + ^B + 6F 

+£2^^ + e^B^^o) _^ g^^(o) ^ (^20) 

= ai (-goo^) , + xQ , (A21) 
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2. Source for the BPT equation 

The source of the BPT equation J2.8t . Sg {r, lo), can be 
expressed in terms of the components of the stress-energy ten- 
sor as follows (see lUsllsUl ') 



S, (r,w) = -V2A]Vr4T(„)(„), 
/2A'AL> 



-^T(n)(fn) e 





2r ^ 



(m)(m) e) 
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where D+ = d/dr^ + iuj, and T'(„)(„) ^ , T(„)(^) ^, r(^)(,-„) ^ 
are the Fourier transform of the following quantities 



T(n)in) e = J dn QSia{d)T'^''n^n^ , (A33) 
r(„)(^) ,= fdn -iS,o{0)T^''n^m, , (A34) 



'-(fn){fh)e = j dn ^2Sm{0)T^"''m^m,^ , (A35) 

where dO. = sin 6d0d4> is the solid angle and the one-forms 
n^, fh^ are defined as 



1 / 2M 
n, ^ -(^1-— ,1,0,0 

= p (0, 0, r, — ir sin^ 

v2 
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Here, -sSem{(^, 4>) ^6 the spin-weighted spherical harmonics 
and we have considered only those with m = because of the 
underlying axisymmetry in the perturbations. 
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